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I. INTRODUCTION 



Fully-developed turbulence is one of the important topics in fluid dynamics. Although it 
has been studied over a century, we have not yet gotten any clear concept of what turbulence 
is. Even if we consider only the simplest case of 'Homogeneous Isotropic Turbulence', we 
do not know what the essential point is. This is mainly because there has been no powerful 
tool to investigate it on either experiment or theory. In the experiments, the data including 
most information is no more than flow lines of trace-particles on the two dimensional sub- 
section, which is far from enough to investigate the fluid in details. Theorists have no 
general method to deal highly nonlinear systems with many degrees of freedom which the 
turbulence belongs to. Even using numerical methods, it is very hard to get the solution 
of fully developed turbulence because of 'fractal property' of it |TJ]. If a system has fractal 
nature, we cannot find any characteristic length. In such an occasion, we must calculate on 
a large system to include many kinds of scale length. This requires us to use large lattice 
sizes (or many number of wave numbers) which needs a long cpu times on computer. Even 
the largest and fastest computers cannot give us the numerical solution of fully developed 
turbulence, although recently some part of characteristic features can be reproducible on 
computer gggg 

In this paper, we propose to give a lattice model which can reproduce some of character- 
istic feature of statistical property of turbulence, ex. scaling properties etc. If we succeed in 
constructing a simple model showing some behaviors which are similar to the fluid, we can 
say that the origin of them is in the common parts. Therefore it will help us to know what 
the essence of turbulence is. 

The paper will be organized as followings. In Sec. [Tl] we will denote on which view point 



we stand and the description of fluid on 'lattice' will be defined. In Sec. Ill, we will define 



a dynamics (time-marching) of our model and show it is equivalent to the Navier-Stokes 
equation. The results and discussion will be presented in Sec. [TV] and Sec. |V| will provide us 
a summary and a conclusion. 



II. MODELING FLUID ON THE LATTICE 
A. Concept 

The concept which we employ in this paper is a vortex tube representation of fluid 
turbulence. This is because it can easily provide us a lattice model to simulate turbulence. 

Recently, in the field of statistical physics, simple lattice models have been shown to be 
good tools to reproduce qualitative features of non-linear systems even if they were essentially 
continuous ones. For example, the diffusion limited aggregation on lattices has turned out 
to be a good model to reproduce fractal dimension of the Laplacian growth systems which 
are essentially continuous ||. K. Kaneko J7| has proposed coupled map lattices in order to 
study chaos in the many degrees of freedom systems. And Y. Oono || has proposed the cell 
dynamical system on lattices and succeeded in getting good results on the phase separation 
process, which has been believed to be described by the Time-Dependent- Gintzburg-Landau 
eq.; a sort of partial differential equation. And one of the author (H.T.) has also proposed 
many discretized models P ]iOl , [TTH and obtained results consistent with the experimental 
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data. These lattice models are more suitable to be treated on computer than conventional 
method like partial differential equations. Following this line we purpose to construct a 
lattice model which can reproduce qualitative features of turbulence in this paper. To 
construct these lattice models, we need some physical quantities whose local dynamics are 
well known. And introducing interactions between them, we can get a complete system to 
represent a real system. For this purpose, the vortex tube representation is very convenient. 
The vortex lines are believed to move with the fluid because of the conservation law of 
circulation. In the following, we will construct the lattice model so as to preserve these 
essential features. The first step of this process is to answer the question; 'How can we map 
the fluid on the lattice?' 

B. The representation of the fluid on the lattice 

To map the fluid on the lattice, we must take special care for some basic values on the 
vector analysis. For example, the divergence of velocity v, divv, must be easy to calculate, 
because incompressibility condition requires divv = 0. And rotv also should be calculated 
easily because of the definition of the vorticity cJ;rotv = uj. Although there is no uniqueness 
on how to map the fluid on the lattice to satisfy these requirements, we employ the following 
description. 

First, we map the velocity field v(r) on each bond on the simple cubic lattice (we will call 
it as the v-lattice), such that x bonds must have the only x component of the velocity, y and 
z bonds also should have the only component along the bond direction. On the other hand, 
the vorticity field u;(r) is assigned on each bond on the dual lattice(we call it ^-lattice). 
Again each bond has only one component. This means, on the v-lattice, the vorticity is 
located on the center of each face and its direction is perpendicular to the face (See Fig. [[]). 
For example, the position vector of the bond between the lattice points (x , yo, z ) and 
(x + 1, yo, Zo) is written as (x + 1/2, y , z ). In this description, the vorticity field ul(r) is 
related to the velocity v(r) by the equation 

Ui(r) = v k {r + ij) + ^(r - ^k) - v k {r - -j) - ^(r + ^k), , (1) 

where i, j, k are unit lattice vector of x, y, z directions, which can be exchanged their order 
cyclically. Therefore uj is circulation rather than the vorticity. However we are not interested 
in the difference between them, because they share the same basic property mentioned in 
the previous subsection, for example folding etc. Independent of whether our cJ(r) is the 
vorticity or the circulation, we get the same description. 

Now, we have everything necessary to construct the map to the lattice model. For 
example, divv(r) is given as 

divv(r)= Yl ± v ( r ±^)> (2) 

i=x,y,z 

where x, y and z are the unit lattice vectors of x, y and z direction respectively. However 
in this description, we have a problem. If we give v(r) first, then we can calculate cJ(r) 
without any problems. But if uj(r) is given first, we can not obtain v(r) uniquely. As will be 
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shown, this is necessary for time evolution of our model. In addition to this, we xrequire the 
incompressibility divv(r) = 0. To get v(r) satisfying this condition seems not to be trivial. 

The basic procedure to calculate v(r) from u(t) is essentially solving a set of linear 
equations given by ( |] ). One may think that it is a similar task to solving the Laplace 
equation. However, this is much more difficult because we cannot use the relaxation method. 
The method which we employ is the eg method which is explained briefly in Appendix A. 
In principle, this method can give us accurate numerical solutions even for large matrices. 
However, if we assume the periodic boundary condition for the lattice, we have another 
problem. Because of the periodicity, one of a set of linear equations mentioned above is 
not independent to the rest equations. Although we cannot prove this point analytically, 
the rank of coefficient matrix of the set of linear equations is always less than the number 
of components by unity. (The total number of the components is equivalent to the total 
number of independent components on the v lattice. This is 3 x L 3 if the linear dimension 
of the lattice is L.) Therefore we have one degree of freedom on ambiguity of values of v(r). 
This ambiguity is used to solve the last problem ; to fix v(r) so as to satisfy divv(r) = 0. 
It looks a difficult task because we have only one degree of ambiguity of the velocity field. 
Fortunately this problem turns out to be solved by requiring a discretized representation 
of divv = at an arbitrary point r . This gives us the coefficient matrix whose dimension 
is equal to its rank. The velocity field v(r) which this matrix gives turns out to satisfy 
divv = and rotv = to. Therefore we succeed in mapping fluid on the lattice. 

III. DYNAMICS OF THE LATTICE MODEL 

In the previous section, we have finished introducing basic elements to the lattice model; 
ctJ(r) and its 'interaction' v(r). Actually, our model corresponds to a set of the vortex tubes 
moving with the velocity induced by the Biot-Savier interaction. In this section, we define 
a concrete process of dynamics. 

In our model, vortex tubes are represented as a sequence of bonds having non-zero 
vorticity. It can be folded and have branches but no ends. This character must be held even 
while the vortex tubes are transferred by the flow of the fluid. The simplest dynamics to 
satisfy these condition is to move the vorticity on some i-bond to the neighboring i-bond 
along j-direction by the amount of u>i(r)vj(r)At and to add some vorticity to neighboring 
j-bonds in order to have no ends. However in stead of that, to move half of the net flow 
forward and the rest half backwards is more favorable in order to recover isotropy (See the 
Appendix B). If we consider i — j makes no change, the time development is defined as; 

u z (r , t + At) = uj z (r , t) + At[±J yz (r =F y, t) ± J xz {r =F x, t) 

TJzx(r - z/2 =p x/2, t) =F J zy (r - z/2 =f y/2, t) 
TJzxijo + 1/2 =F x/2, t) =F J zy (r + z/2 t y/2, t)]/2 

(3) 

where Jy(r, t) = Vi(r,t)ujj(r,t) and Vi(r) is the i-th component of the velocity on the u- 
lattice and is defined as the average over neighboring bonds in the v-lattice. For example, 
v on z-bond is defined as 
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v(r„) = {[v x {r + f/2) + v x {t - y/2)]/2, [v y {r + k/2) + v y {r - x/2)]/2, 0) 

7/ and z components have similar time developments. At is the time interval and chosen so as 
to max r | Jij(r, t)At \— a. This makes the change of the strength of u)(r) less than a where 
a is a small positive constant. In addition to the above, we add the renormalization of v(r) 
and u(t) at each step. The above procedure does not guarantee the conservation of energy. 
Although the increasing of energy at each step is small, it may cause some influence for 
long time behavior. For safety, we renormalize v(r) and w(r) just after the vortex dynamics 
given by ( | ). 

Next we introduce the viscosity. The viscosity is introduced as a diffusion process of the 
vorticity. And it is easily represented as; 

to z (r ,t + At) = {1 - vAt)u z {r ,t) + ]T u) z (r ±i,t)vAt/6 (4) 

i=x,y,z 

where v is the viscosity, y and z components follow similar equations. This is the conven- 
tional and the simplest discretization of the diffusion equation. 

In addition to the above, if we would like to maintain steady state, we should introduce 
an external force. The time development of u(r) due to the force is very simple, 

uj z (r 0l t + At) = uj z (r , t) + rotf (r )At (5) 

where f is the external force. 

The whole process of our algorithm is summarized as follows; 

1. Calculating £ r v(r) 2 /2 

2. Deciding At 

3. Time development corresponding to ( ^ ) 

4. Renormalizing v(r),cD(r) so as to keep X)r v ( r ) 2 /2 constant. 

5. Time development corresponding to the diffusion process. (See ( [| )) 

6. Time development corresponding to the external force. (See ( [5] )) 

7. Return to process 1. 



IV. THE NUMERICAL RESULTS 

In this section, we show some numerical results. In the numerical simulation, we use the 
lattice with L = 24 and set a=0.1. External force f(r) is given as 

f x (r) = 4 sm(k(x + -)) cos(ky) cos(kz) 

1 1 

+ 2(cos(/c(x + -)) cos(ky) + cos(ky) cos(kz) + sin(A;(x + -)) sin(A;2;)) 

f y (r) = —4 cos(kx) sin(k(y + -)) cos(kz) 

+ 2(sin(k(y + -)) sin(A;x) + cos(&;,2) cos(k(y + -)) + cos(&;,2) cos(A;x)) 
f z (r) = 2(cos(k(z + -)) cos(A;x) + cos(A;x) cos(ky) + sin(%) sin(k(z + -))) 
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where k = 2tt/L. Actually, these correspond to the Taylor-Green vortex accompanied with 
the components of wave number 2\^2n/L. This choice is not unique but, this is the case 
which contains the fewest number of kinds of wave number among those which recover 
isotropy in the numerical simulations. 

Remaining parameter is the viscosity v. We investigate four values of v: 5, 8, 10, 20 x 10 -3 . 
If we define microscopic Reynolds number R\ = y / 10/3(v 2 /2)/z/y / (u 2 /2), R\ = 5 ~ 20 ((• • •) 
means average over the whole lattice point and time). However we are not sure whether 
absolute values have meanings or not, because in our model the unit cubic has already been 
coarse grained. 

First we show the time development of the total I] r v(r) 2 and J2r^{ T ) 2 m Fig 0- These 
behaviors look like those of fully developed turbulece gotten with the conventional direct 
integration (For example, Ref. |12[| ). Every results shown in below are averaged over nine 
snap shots apart from each other by the same intervals after the state reaches steady state 
(Actually speaking, 600 < t < 3000). One can see total square velocities are not influenced 
by the change of v. On the other hand, the total square vorticities are influenced by it 
strongly. 

Before going forward, we show snapshot configurations at a instance. Fig. |] shows a 
velocity field snapshot. Nothing abnormal is observed in this figure. In Fig. [5], the bonds 
on which strong vorticities exist are plotted. The number of plotted bonds are 527 which is 
much less than total number of bonds : 41472. In spite of small ratio of number of plotted 
bonds, we can see large number of connected bonds. This means our velocity field surely 
succeeds in presenting coherent structure of vorticities which is observed in direct integration 

For each calculations with four different z/, we need less than nine hours on CRAY-XMP 
(3000 steps each). This is surprisingly short, if obtained velocity fields really correspond to 
fully developed turbulence. 

To check this point, three dimensional energy spectrums are shown in Fig|| One must 
take care of that drastic decreasing of energy spectrums on high wave number region is due 
to not only the energy dissipation but also the decreasing of number of mode having wave 
number of k. So we should consider only | k |< 7T. In order to see whether they exhibit 
inertial range, that is, k~ 5 ^ 3 law, the averaged spectrum is shown in Fig. [7]. For u — 10 
and 20 x 10~ 3 , we can see inertial ranges. However, for other two with smaller viscosity, 
humps are seen in higher wave number region. This means, for later two cases, dynamics 
in small scale are damaged and they cannot be regarded as a good representation as highly 
developed turbulence. 

One may think these spectrum do not have any dissipation ranges, because decreasing 
of energy spectrum may be caused by only insufficient number of modes as pointed out in 
above. For these four cases, we also calculate one dimensional energy spectrum E z =<\ 
J v z (x,y, z)exp(—ikx)dx \ 2 > y , z , where < • • • > y<z stands for the average over y, z. We can 
see some scaling region obeying Kolmogorov's scaling region and short inertial subrange. 
Here we can see clear dissipation range. This time decreasing of energy spectrum at the 
higher wave number is surely caused by the energy dissipation because shortness of number 
of modes does not occur in one dimension. 

Therefore, we think we can think cases with two larger viscosities can be investigated 
as representing fully developed turbulence, but not for other two cases. However, we show 
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results for all four cases for comparison. 

In Figs|| the probability distribution functions (PDF) of v x are shown. They look like 
almost Gaussian which is also observed in experiments and numerical calculation JI. In 
order to see the dependence upon z/, semi- logarithmic plot of them is shown in Fig.[lC. 
Those having smaller v are closer to Gaussian than others with larger v. This means even 
if higher wave number components are damaged, lower component can properly behave and 
become closer to fully developed turbulence. This also reveals that our model works well to 
produce developed turbulence. 

Next several PDFs which are close to exponential-like distributions are shown. Figs.|TT| 
to [L4| show PDF of lu x , dv x /dx,dv x /dy and v x having only higher wave number components 
k > 27r/3. All of them are normalized so as to have standard deviations of unity. Roughly 
speaking all of them look like exponential distribution. They agree with the results obtained 
before p|,|^,|^,[5| Jl^Jl^Jl^JT6[| well. Especially, asymmetry of PDF of dv x /dx is properly repro- 
duced. It is extremely surprising because these PDFs reflect the property in the high wave 
number region which is damaged for lower viscosities. The fact means that our model can 
catch the essential feature of turbulence which is hardly destroyed even by the insufficient 
resolution. It enables us to get some aspects on what the essence of turbulence is fL7| . 

These results clearly demonstrate our model works so well. It is clear that our model can 
work very well to represent fully developed turbulence although it requires small memory- 
sizes and short cpu-time. 



V. CONCLUSION AND SUMMARY 

From the above results, our model works remarkably well. For example, since our model 
has succeeded in reproducing some characteristic features of fully developed turbulence; 
the Kolmogorov's scaling of the energy spectrum, the development of energy and enstrophy, 
coherent structure of vortex tubes, the qualitative feature of the distribution functions of the 
velocity and its differential we can conclude that we have achieved our purpose; to construct 
the simple model which can reproduce some qualitative features of turbulence. We think 
it proves that the lattice model is also valid in investigating the turbulence. In contrast to 
the conventional direct integration of Navier-Stokes equation, our modeling is very tough to 
reproduce the qualitative features of turbulence. Our model is not damaged even with the 
crude procedures like artificial rescaling in order to conserve energy and rough discretization 
of the field, which cause severe catastrophe to the direct integration. 

Moreover we can expect many things for our model. Since our model is so simple, we can 
find some essential reason about the appearance of the non-Gaussian PDFs, which will be 
reported elsewhere |17[]. And even some rigorous calculation like estimating upper bounds 



of some value may be possible, because many such calculation have been done on the lattice 
systems. In addition to this, our model is connected with other lattice models in statistical 
physics. This will allow us to investigate several phenomena like chemical reaction and phase 
ordering process under the existence of turbulent flow. 
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APPENDIX A: THE CG METHOD 

In this appendix, we give a brief explanation of the eg (conjugate gradient) method. For 



details, see references |T8|,|T9 



Assume we have a set of linear equations. In matrix form, Ak = b, where A is the 
coefficient matrix, x is a vector consisting of variables and b is the constant vector having 
the same length as that of x. The following procedures can be applied to the case that A is 
not symmetric neither positive definite. 

The eg method is one of those which are suitable to search extreme values. If one would 
like to search the extreme values of the function /(r), one can follow along the flow V/(r) 
iteratively. These are applicable to find the solution A~k = b, if we define /(r) as r 2 with 
r = A~k — b. To accelerate this process, in each iteration of the eg method, the direction 
of flow is tuned so as to be orthogonal to the flow directions in all the previous iterations. 
Therefore the eg method should converge within the n iterations, where n is the number of 
components; that is, the total number of dimension in the variable space. Skipping detailed 
points, the actual algorithm of eg method is 

1. For given fe-th approximate solution x^ and fc-th residue r^, prepare a store array 
Pfc. For k=l, pi = i"i = b — Axi 

2. calculate a k - {ApkApk) - {Apk ,A Pk ) 

3. Calculate the k + 1-th approximate solution x^+i = x^ + a k p k 

4. Calculate the k + 1-th residue r^+i = b — Ax^i = — a k Ap k 

Pfllpulfl+P R, - {A-y k ,AA T r k+1 ) _ (A T r k+1 ,A T r k+1 ) 

{j. calculate p k - {Apk ,Ap k ) ~ (^r k ,A^r k ) 
6. Calculate p k +i = A T r k+1 + (3 k p k 

For the convenience of numerical simulation, one can do as followings using store vectors 
g, q. For k — 1, i"i = b — Axi, p = 0, c = 1, and iterate like 

g fc = A T r k 

c k = (g fc , gfc) 

Pk = c k /c k -i 

Pk = gfc + PkPk-i 

q fe = Ap k 

4 = (q*, q*) 

«fc = c k /d k 

X fe+ l = X fe + OfcPfc 

r fc+ i = r k - a k q_ k 



S 



In our case, x^ consists of each component of velocity Vi(r). Therefore if the linear 
dimension of the lattice is L, the total number of the component n is 3 x L 3 . A is coefficient 
matrix of ( [1] ) and ( )• b is the vector containing all the components of u(r). 



APPENDIX B: RELATION WITH THE VORTICITY EQUATION 



One may think that our model lacks mathematical justification and therefore cannot be 
trust. In principle, we do not think our model needs any justifications, because anyway our 
model is purposed only to reproduce qualitative results at the moment. If one has such a 
purpose, one dose not have to make his model being related with the real system exactly. 
This attitude is justified with the hope of existence of 'universality class'. The universality 
class is a set of model or parameters which can give the same results on some property. 
One of the famous universality class in statistical physics is 'Ising universality class' in the 
magnetic phase transition. Ising model is a set of magnetic moments which takes only 
values of ±1, therefore far from the true magnetic materials. However, surprisingly, it can 
reproduce almost all basic nature of 'Critical Phenomena'. Although no one cannot prove a 
true magnetic material is exactly same as the Ising model, one can get many deep insights 
about the 'true' phase transition accompanied with the critical phenomena. And it has 
helped many people who would like to know what the critical phenomena are. 

We stand at the same point on the turbulence. Therefore if our model can reproduce 
some qualitative feature of turbulence, it is useful enough to use in order to investigate 
turbulence. In this meaning, we do not need any special justifications other than the results 
itself. However, we can actually shown some relationship with the true fluid equation as 
follows. 

In order to expand ( |3] ) in space, first we expand v(r) and uj(t) itself. 

00 1 

F,(r + 5± + 5y + Sz) = £ -r(&— + 5y— + 6z—) n Fi 

£t q n\ Ox Oy Oz 


= exp(5x— + 5y— + 5z—)Fi 
ox Oy oz 



where Fi is either v$ or tUj. Considering the definition of v,we get from eq.(3), 
Ui(t + At) = Ui(t) 

5j 8% 

+ At[cosh(y— ) sinh(5fc— )(-v k u i + 2 cosh(y — )viU k ) 



: yi )sinh( ^ )( -^ ! + 2cosh( 2a ; 



where i,j, k is x,y,z respectively and change order cyclical and the partial differentiation 
with tilde is applied to only the velocity not to the vorticity. And taking into account to 
the second order of Si, Sj, Sk and the first order of 5t, we can get 

duAt) r .d, . d 

gt = Qj — iViUJj - VjUJi) + bk—{ViUJ k - v k 0Ji) 
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Using divv = divu; = and setting 5x = 5y = 5z = 1, this reduces to 

^ + ( v • V)v = {Co ■ V)v. 

This is the vorticity equation without viscosity term itself. And it does not have 'any second 
order' terms. This means we can introduce viscosity term using eq.( |] ) correctly. Therefore 
our model reduces to the dynamics of vortex by taking the continuum limit. 

Here we should mention our dynamics differ from the conventional ones. In the conven- 
tional concepts, vortex tube should move with fluid in nonviscostic case. In our dynamics, as 
show in Fig. |3| the center of vortex tube does not move. In stead of that, the width |2D[ of the 
vortex tube will increase rapidly. We see that ours are better than conventional ones. This 
is because the conventional ones give the second order term which vanish in our dynamics. 
The second order terms conflict with the true viscosity term which is also of the second 
order. We guess this is the reason why the conventional vortex method cannot either get 
good results nor consider viscosity correctly. Therefore if one follows our dynamics, vortex 
method can give us physical results. 
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FIGURES 

FIG. 1. The relation of a vortex tube (the doubled vector) and velocity vectors (the solid 
vectors), u is equal to v\ y — V2 y — V2 Z + v\ z . 

FIG. 2. The actual changes of the distribution of "vortex tubes" on the w-lattice during the 
dynamical process. If the part of the vortex tube has the vorticity uo and the velocities on that part 
are v x , v y (left), after the dynamics it has branches of vortex tubes shown as narrow arrows(right). 
The direction of the arrows shows the sign of the value of the vorticity. And absolute values of 
the branches have two values depending on whether it was induced by v x or v y . One must notice 
the vorticities on the part which has existed before dynamics does not change at all. This shows 
in our dynamics the canter of the vortex tube does not move and the vortex tube itself collapses 
automatically. In Appendix. B we can see this is better than the conventional picture. 

FIG. 3. Time development of v(r) 2 (thin lines) and w(r) 2 (thick lines) (a) v = 5 X 10~ 3 (b) 
v = 8 X KT 3 (c) v = 10 x 1(T 3 (d) v = 20 x 10~ 3 

FIG. 4. A velocity snap shot taken from v = 8 X 10~ 3 at t ~ 1500. Two dimensional cross 
section at z = 4. 

FIG. 5. Coherent structure of vortex tubes at the same instance as that in Fig. |j. Bonds with 
strong vorticities are plotted. Summation over them is 14 percent of total square vorticities. 

FIG. 6. Three dimensional energy spectrums. (a) to (d) correspond to those in Fig. ||[ Solid 
lines indicate A; -5 / 3 . 

FIG. 7. Averaged spectrum corresponding to Figs. |6| Whole range is divided into 50 parts and 
spectrum is averaged within each band. 

FIG. 8. The scaled energy spectrum, = (e/z^ 3 ) 1 / 4 7r/12 is proportional to Kolmogorov wave 
number, e = v(uj 2 /2) is dissipation energy. Solid line indicates 5/3 law. The symbols represent 
the cases with different values of the viscosity rate v\ x : 5 x 10 -3 , □ : 8 x 10~ 3 , O : 10 x 10 -3 , 
+ : 20 x 10~ 3 . Error bars are smaller than the size of character. 

FIG. 9. PDF of v x . (a) to (d) correspond to those of Fig. [||. The lines without chracter 
represent Gaussian distribution. 

FIG. 10. Semi- logarithmic plot of Fig. ||. x,D,0 and + correspond to (a),(b),(c) and (d), 
respectively. 



FIG. 11. PDF of uj x . Characters are the same as those of Fig. 10 
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FIG. 12. PDF of dv x /dx. Characters are the same as those of Fig,|l~0| 



FIG. 13. PDF of dv x /dy. Characters are the same as those of Fig.lC 



FIG. 14. PDF of v x . Their lower wave number components k < 2tt/3 are removed. Characters 



are the same as those of Fig. 10 
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